
*** This do file includes regressions based on the National Risk and Vulnerability Assessment (NRVA) survey for the paper "Foreign Interventions and Community Cohesion in Times of Conflict"	
*** Author: Sarah Langlotz		
*** Journal: Journal of Development Economics
													

***********************************************************************************************
**************************************Read me **************************************
***********************************************************************************************
*** FOLDER STRUCTURE
/* After downloading the data from the Harvard Dataverse, create the following three folders:
	- "processed"
	- "graphs"
	- "tables"
*/



***	DATA ACCESS
*	The raw NRVA household survey data are subject to access restrictions and cannot be shared publicly. 
*	Please contact the author for information on how to request access.



***	SOFTWARE VERSION
*	Stata/MP 19.0 for Windows (64-bit x86-64), Revision 21 May 2025

*	Please note that rdrobust produces different bandwidths depending on the version. This paper used *!version 7.3.0  09Jun2017  
net install st0366_1, from(http://www.stata-journal.com/software/sj17-2) replace



***	TABLES AND FIGURES
* 	This dofile creates Figures 2, 3, 4, 5
* 	This dofiles created Tables 1, 2, 7


***********************************************************************************
***	Set your directory

local DIR = ""
cap cd "`DIR'"

***********************************************************************************

drop _all
clear matrix
clear mata
set maxvar 10000
set matsize 10000

set more off

clear
use "processed\finaldata_hhlevel_NRVA_external.dta"
	
	
	global rddcontrols anyshock l1lnaid_wb l1vhi log_l1nightlight loan_i l1lnallsigacts
	global rddcontrols_short l1lnaid_wb l1vhi log_l1nightlight  l1lnallsigacts
	global short l1lnaid_wb l1vhi log_l1nightlight loan_i
	

********************************************************************************
* FIGURE 2 - GRD: Balancing tests
********************************************************************************

		cap drop indivbw
		gen indivbw=.
				
		cap est drop *
		
		preserve
		*balancing in coef plots for standardized outcomes
		foreach Y of varlist anyshock shock amount_wheat_flour lnhhexp_2005 insecurity copehelp copesocial ///
		sumassets_2003 economicimprove2003 copemilstate ///
		share_native1 share_native2 share_native3 kuchi2003 hhmember2003 {
		qui rdbwselect `Y' fdist_disoldb dum_segment200km* if year==2003 &  province!="bamyan" & province!="kunduz" & province!="balkh" & province!="paktya" & province!="kabul", c(0) p(1) kernel(uniform) vce(cluster id) 
		replace indivbw=e(h_mserd)
		
		qui reghdfe `Y' treat_old fdist_disoldb [pweight=hhweight] if year==2003 & dist_disoldb<indivbw &  province!="bamyan" & province!="kunduz" & 		province!="balkh" & province!="paktya" & province!="kabul", absorb(segment200km) cluster(districtcode)
		sum `Y' if  e(sample) 
		gen `Y'_sd=(`Y'-r(mean))/r(sd)
		sum  `Y'_sd if  e(sample)

		}

				
		* Analysis: all coefficients
		mat B=J(15,2,0)
		local no = 1

		estimates drop _all	

		foreach Y in anyshock_sd insecurity_sd shock_sd copehelp_sd copesocial_sd  ///
		 copemilstate_sd amount_wheat_flour_sd lnhhexp_2005_sd sumassets_2003_sd economicimprove2003_sd ///
		share_native1_sd share_native2_sd share_native3_sd kuchi2003_sd hhmember2003_sd {
		rdbwselect `Y' fdist_disoldb dum_segment200km* if year==2003 &  province!="bamyan" & province!="kunduz" & province!="balkh" & province!="paktya" & province!="kabul", c(0) p(1) kernel(uniform) vce(cluster id) 
		replace indivbw=e(h_mserd)
				
		 eststo est`Y': reghdfe `Y' treat_old fdist_disoldb [pweight=hhweight] if year==2003 & dist_disoldb<indivbw &  province!="bamyan" & province!="kunduz" & province!="balkh" & province!="paktya" & province!="kabul", absorb(segment200km) cluster(districtcode)
		estadd local indvw=round(indivbw, .1)
 scalar b`no' = _b[treat_old]
		 scalar s`no' = _se[treat_old]
		 matrix B[`no', 1] =  b`no'
		 matrix B[`no', 2] =  s`no'
		 local no = `no'+1
		 
		}
		
									
		* Draw and export
		coefplot matrix(B[,1]), se(B[,2]) ///
    levels(99 95) ///
    ciopts(lcolor(gs0) fcolor(gs10)) /// 
    xsize(13) ysize(16) ///
    xline(0, lwidth(thin)) scale(0.75) ///
    coeflabels(r1="Any shock" r2="Insecurity shock" r3="Climate shock" r4="Comm. help" r5="Comm. help+loan" r6="Cope state/military" ///
               r7="Wheat consumption" r8="Food expenditures" r9="Sum of assets" r10="Economic improve" ///
               r11="Language: Dari" r12="Language: Pashto" r13="Language: Uzbeki" r14="Kuchi" r15="Number members") ///
    xtitle("Effect on standardized outcome", color(black) margin(medium)) ///
    mcolor(gs0) msize(medium) /// 
    graphregion(color(white)) ///
    xlabel(-1.5(0.5)1, nogrid labgap(1)) ///
    headings(r1="{bf:Panel A: Shocks, coping}" r6="{bf:Panel B: Income}" r11 = "{bf:Panel C: Ethnicity, size}")


		graph export "graphs\coefplot_balancinghh.pdf", replace				
		graph export "graphs\coefplot_balancinghh.png", replace				
		
	restore	
		
		
		*collapse at district-level 
		cap drop control1996
		encode groupcontrol_taliban1996, gen(control1996)
		preserve 
		keep if year==2003
					
						collapse (mean) treat_old fdist_disoldb dist_disoldb l1lnbest_est l1lnallsigacts  l1aid_afg l1aid_wb ruggedness_1000 suitability_wheat ///
						revenue_tons travel_mean rural_share control1996 pashtun_ethno pashtun_greg no_ethnic_greg   pop log_l1nightlight activebase largebase_inbr suitability_opium (last) wave province segment200 , by(districtcode)
						gen lnpop=log(pop+1)
						gen indivbw=.
						
						cap est drop *
							
		
		
		*balancing in coef plots for standardized outcomes

		foreach X of varlist l1lnbest_est l1lnallsigacts activebase l1aid_afg l1aid_wb log_l1nightlight  ruggedness_1000 suitability_wheat suitability_opium lnpop travel_mean control1996 pashtun_ethno pashtun_greg no_ethnic_greg {
		rdbwselect `X' fdist_disoldb segment200kmdum* if province!="bamyan" & province!="kunduz" & province!="balkh" & province!="paktya" & province!="kabul", c(0) p(1) kernel(uniform) 
		replace indivbw=e(h_mserd)
		qui reghdfe `X' treat_old fdist_disoldb if dist_disoldb<indivbw &  province!="bamyan" & province!="kunduz" & province!="balkh" & province!="paktya" & province!="kabul", absorb(segment200) vce(robust)
		sum `X' if  e(sample) 
		gen `X'_sd=(`X'-r(mean))/r(sd)
		sum  `X'_sd if  e(sample)

		}

				
		* Analysis: all coefficients
		mat C=J(15,2,0)
		local no = 1

		estimates drop _all	

		foreach Y in l1lnbest_est_sd l1lnallsigacts_sd activebase_sd l1aid_wb_sd l1aid_afg_sd  log_l1nightlight_sd  ruggedness_1000_sd suitability_wheat_sd suitability_opium_sd lnpop_sd travel_mean_sd control1996_sd pashtun_ethno_sd no_ethnic_greg_sd {
		rdbwselect `Y' fdist_disoldb segment200kmdum* if province!="bamyan" & province!="kunduz" & province!="balkh" & province!="paktya" & province!="kabul", c(0) p(1) kernel(uniform) 
		replace indivbw=e(h_mserd)
		eststo est`Y': reghdfe `Y' treat_old fdist_disoldb if dist_disoldb<indivbw &  province!="bamyan" & province!="kunduz" & province!="balkh" & province!="paktya" & province!="kabul", absorb(segment200) vce(robust)
		estadd local indvw=round(indivbw, .1)
		 scalar b`no' = _b[treat_old]
		 scalar s`no' = _se[treat_old]
		 matrix C[`no', 1] =  b`no'
		 matrix C[`no', 2] =  s`no'
		 local no = `no'+1
		}
		
		* Draw and export
		coefplot matrix(C[,1]), se(C[,2]) ///
    levels(99 95) ///
    ciopts(lcolor(gs0) fcolor(gs10)) /// 
    xsize(13) ysize(16) ///
    xline(0, lwidth(thin)) scale(0.75) ///
    coeflabels(r1="UCDP BRD" r2="SIGACTS events" r3="Active bases" r4="WB aid" r5="Total aid" ///
    r6="Log night light" r7="VHI" r8="Ruggedness" r9="Wheat Suitability" r10="Opium Suitability" ///
    r11="Log Population" r12="Travel time" r13="Territory control" r14="Pashtuns" r15="Number ethnic groups") ///
    xtitle("Effect on standardized outcome", color(black) margin(medium)) ///
    mcolor(gs0) msize(medium) ///
    graphregion(color(white)) ///
    xlabel(-5(1)5, nogrid labgap(1)) ///
    headings(r1="{bf:Panel A: Conflict, aid}" r6 = "{bf:Panel B: Development, geography}" r11 = "{bf:Panel C: Socio-demographics}") 


		graph export "graphs\coefplot_balancingdistrict.pdf", replace				
		graph export "graphs\coefplot_balancingdistrict.png", replace				
		
restore						
		
		
		
********************************************************************************
* FIGURE 3 - Panel: Reconstruction, formal institutions, and development, 2005--2010
********************************************************************************

	
	local mech_outcomes /// 
	copemilstate loan agrincome eradication /// 
	f_enrolled m_enrolled sh_kids_enroll sh_f_enroll sh_f_m_enroll /// 
	amount_wheat_flour lnhhexp_totp2005 dietarydiversity foodsecure sumassets /// 
	cashforwork lnaid_wb log_nightlight economicimprove 
	
	
	foreach Y of local mech_outcomes {

    * check variable exists
    capture confirm variable `Y'
    if _rc {
        di as error "Variable `Y' not found – skipped"
        continue
    }

    * compute mean and sd on analysis sample
    quietly summarize `Y' if year>2003 & year<2011

    * skip if variable has zero variation
    if r(sd)==0 {
        di as error "Variable `Y' has zero SD – skipped"
        continue
    }

    * create standardized version
    gen z_`Y' = (`Y' - r(mean)) / r(sd) if year>2003 & year<2011
    label var z_`Y' "`: var label `Y'' (standardized)"
}
 

cap grstyle init
cap grstyle color background white
cap mkdir figures


local outcomes ///
    copemilstate loan agrincome eradication ///
    f_enrolled m_enrolled sh_kids_enroll sh_f_enroll sh_f_m_enroll ///
    amount_wheat_flour lnhhexp_totp2005 dietarydiversity foodsecure sumassets ///
    cashforwork lnaid_wb log_nightlight economicimprove


local dist_outcomes "log_nightlight lnaid_wb"


local samp "year>2003 & year<2011"


quietly summarize no_securityprt if `samp'
local sd_prt = r(sd)
if missing(`sd_prt') | `sd_prt'==0 {
    di as error "SD of no_securityprt is missing/zero in the estimation sample. Cannot SD-scale."
    exit 459
}
di as txt "SD(no_securityprt) in sample = " as res `sd_prt'


local fe_hh "id i.year"
local fe_d  "districtcode i.year"

*==============================================================================*
* A) CONTEMPORANEOUS
*==============================================================================*
local K : word count `outcomes'
matrix A0 = J(`K',2,.)
matrix B0 = J(`K',2,.)
matrix C0 = J(`K',2,.)
local coeflabs0 ""

local no = 1
foreach y of local outcomes {

    * Y-axis label
    local lab : variable label `y'
    if "`lab'"=="" local lab "`y'"
    local coeflabs0 `"`coeflabs0' r`no'="`lab'""'

    *==========================
    * DISTRICT-LEVEL outcome
    *==========================
    if strpos(" `dist_outcomes' ", " `y' ") {

        preserve
            keep if `samp'

            * Collapse to district-year with needed regressors + controls
            collapse (mean) ///
                `y' activebase largebase_inbr no_securityprt ///
                , by(districtcode year)

            * Standardize outcome on district-year sample (create z_`y')
            quietly summarize `y'
            gen z_`y' = (`y' - r(mean))/r(sd)

            * Any base
            quietly reghdfe z_`y' activebase, ///
                absorb(`fe_d') cluster(districtcode)
            matrix A0[`no',1] = _b[activebase]
            matrix A0[`no',2] = _se[activebase]

            * Large base
            quietly reghdfe z_`y' largebase_inbr, ///
                absorb(`fe_d') cluster(districtcode)
            matrix B0[`no',1] = _b[largebase_inbr]
            matrix B0[`no',2] = _se[largebase_inbr]

            * PRT security projects (scale to 1 SD)
            quietly reghdfe z_`y' no_securityprt, ///
                absorb(`fe_d') cluster(districtcode)
            matrix C0[`no',1] = _b[no_securityprt]  * `sd_prt'
            matrix C0[`no',2] = _se[no_securityprt] * `sd_prt'

        restore
    }

    *==========================
    * HOUSEHOLD-LEVEL outcomes
    *==========================
    else {
				
        cap confirm variable z_`y'
        if _rc {
            di as error "Missing z_`y' -> standardize this outcome first."
            exit 198
        }

        quietly reghdfe z_`y' activebase ///
            if `samp' [pweight=hhweight], absorb(`fe_hh') cluster(districtcode)
        matrix A0[`no',1] = _b[activebase]
        matrix A0[`no',2] = _se[activebase]

        quietly reghdfe z_`y' largebase_inbr ///
            if `samp' [pweight=hhweight], absorb(`fe_hh') cluster(districtcode)
        matrix B0[`no',1] = _b[largebase_inbr]
        matrix B0[`no',2] = _se[largebase_inbr]

        quietly reghdfe z_`y' no_securityprt ///
            if `samp' [pweight=hhweight], absorb(`fe_hh') cluster(districtcode)
        matrix C0[`no',1] = _b[no_securityprt]  * `sd_prt'
        matrix C0[`no',2] = _se[no_securityprt] * `sd_prt'
    }

    local no = `no' + 1
	
}

coefplot ///
    (matrix(A0[,1]), se(A0[,2]) label("Any base")          pstyle(p2) msymbol(O) mcolor(maroon) offset(0.35) ciopts(color(maroon maroon)) level(99 95)) ///
    (matrix(B0[,1]), se(B0[,2]) label("Large base")        pstyle(p5) msymbol(D) mcolor(navy)   offset(0)    ciopts(color(navy navy))     level(99 95)) ///
    (matrix(C0[,1]), se(C0[,2]) label("PRT sec. projects") pstyle(p8) msymbol(T) mcolor(teal)   offset(-0.35) ciopts(color(teal teal))   level(99 95)) ///
    , ///
    xtitle("Effect on standardized outcome", margin(medium)) ///
    xline(0, lcolor(cranberry) lwidth(thin) lpattern(solid)) ///
    coeflabels(`coeflabs0') ///
    headings( ///
        r1  = "{bf: Employment & coping}" ///
        r5  = "{bf: Education}" ///
        r10 = "{bf: Living standards}" ///
        r15 = "{bf: Aid & improvement}" ///
    ) ///
    legend(position(7) rows(1) region(lstyle(black))) ///
    xsize(12) ysize(12) ///
    scale(0.55) ///
    name(mech_current, replace)

graph save   "graphs\mechanisms_current_coefplot.gph", replace
graph export "graphs\mechanisms_current_coefplot.pdf", replace


*==============================================================================*
* B) LAGGED (t-1)
*==============================================================================*
matrix A1 = J(`K',2,.)
matrix B1 = J(`K',2,.)
matrix C1 = J(`K',2,.)
local coeflabs1 ""

local no = 1
foreach y of local outcomes {

    * Y-axis label
    local lab : variable label `y'
    if "`lab'"=="" local lab "`y'"
    local coeflabs1 `"`coeflabs1' r`no'="`lab'""'

    *==========================
    * DISTRICT-LEVEL outcome
    *==========================
    if strpos(" `dist_outcomes' ", " `y' ") {

        preserve
            keep if `samp'

            * Collapse to district-year with lagged RHS vars + controls
            collapse (mean) ///
                `y' l1activebase l1largebase_inbr l1no_securityprt ///
                , by(districtcode year)

            * Standardize outcome on district-year sample
            quietly summarize `y'
            gen z_`y' = (`y' - r(mean))/r(sd)

            * Any base (lagged)
            quietly reghdfe z_`y' l1activebase, ///
                absorb(`fe_d') cluster(districtcode)
            matrix A1[`no',1] = _b[l1activebase]
            matrix A1[`no',2] = _se[l1activebase]

            * Large base (lagged)
            quietly reghdfe z_`y' l1largebase_inbr, ///
                absorb(`fe_d') cluster(districtcode)
            matrix B1[`no',1] = _b[l1largebase_inbr]
            matrix B1[`no',2] = _se[l1largebase_inbr]

            * PRT security projects (lagged), scale to 1 SD of contemporaneous no_securityprt
            quietly reghdfe z_`y' l1no_securityprt, ///
                absorb(`fe_d') cluster(districtcode)
            matrix C1[`no',1] = _b[l1no_securityprt]  * `sd_prt'
            matrix C1[`no',2] = _se[l1no_securityprt] * `sd_prt'

        restore
    }

    *==========================
    * HOUSEHOLD-LEVEL outcomes
    *==========================
    else {
		
        cap confirm variable z_`y'
        if _rc {
            di as error "Missing z_`y' -> standardize this outcome first."
            exit 198
        }

        quietly reghdfe z_`y' l1activebase  ///
            if `samp' [pweight=hhweight], absorb(`fe_hh') cluster(districtcode)
        matrix A1[`no',1] = _b[l1activebase]
        matrix A1[`no',2] = _se[l1activebase]

        quietly reghdfe z_`y' l1largebase_inbr  ///
            if `samp' [pweight=hhweight], absorb(`fe_hh') cluster(districtcode)
        matrix B1[`no',1] = _b[l1largebase_inbr]
        matrix B1[`no',2] = _se[l1largebase_inbr]

        quietly reghdfe z_`y' l1no_securityprt ///
            if `samp' [pweight=hhweight], absorb(`fe_hh') cluster(districtcode)
        matrix C1[`no',1] = _b[l1no_securityprt]  * `sd_prt'
        matrix C1[`no',2] = _se[l1no_securityprt] * `sd_prt'
    }

    local no = `no' + 1
	
	
}

coefplot ///
    (matrix(A1[,1]), se(A1[,2]) label("Any base")          pstyle(p2) msymbol(O) mcolor(maroon) offset(0.35)  ciopts(color(maroon maroon)) level(99 95)) ///
    (matrix(B1[,1]), se(B1[,2]) label("Large base")        pstyle(p5) msymbol(D) mcolor(navy)   offset(0)     ciopts(color(navy navy))     level(99 95)) ///
    (matrix(C1[,1]), se(C1[,2]) label("PRT sec. projects") pstyle(p8) msymbol(T) mcolor(teal)   offset(-0.35) ciopts(color(teal teal))     level(99 95)) ///
    , ///
    xtitle("Effect on standardized outcome", margin(medium)) ///
    xline(0, lcolor(cranberry) lwidth(thin) lpattern(solid)) ///
    coeflabels(`coeflabs1') ///
    headings( ///
        r1  = "{bf: Employment & coping}" ///
        r5  = "{bf: Education}" ///
        r10 = "{bf: Living standards}" ///
        r15 = "{bf: Aid & improvement}" ///
    ) ///
    legend(position(7) rows(1) region(lstyle(black))) ///
    xsize(12) ysize(12) ///
    scale(0.55) ///
    name(mech_lagged, replace)

graph save   "graphs\mechanisms_lagged_coefplot.gph", replace
graph export "graphs\mechanisms_lagged_coefplot.pdf", replace



********************************************************************************
* FIGURE 4 - GRD: Heterogeneity in aid effectiveness, 45km bandwidth, 2005
********************************************************************************	

	global rddaidcontrol anyshock l1vhi loan_i
	
	gen aid=.
	gen treataid=.
						
	cap gen mainbw=45
			

qui reg log_nightlight c.l1lnaid_wb##i.treat_old ///
    $rddaidcontrol fdist_viloldb seg_viloldb2dum* ///
    if year==2005 & dist_viloldb<45 & province!="bamyan" & province!="paktya" & province!="kabul", ///
    cluster(districtcode)
sum l1lnaid_wb if e(sample)
cap drop l1lnaid_wb_sd
gen double l1lnaid_wb_sd = l1lnaid_wb * r(sd)

qui reg log_nightlight c.l1lnaid_afg##i.treat_old ///
    $rddaidcontrol fdist_viloldb seg_viloldb2dum* ///
    if year==2005 & dist_viloldb<45 & province!="bamyan" & province!="paktya" & province!="kabul", ///
    cluster(districtcode)
sum l1lnaid_afg if e(sample)
cap drop l1lnaid_afg_sd
gen double l1lnaid_afg_sd = l1lnaid_afg * r(sd)

foreach Y of varlist amount_wheat_flour lnhhexp_totp2005 dietarydiversity foodsecure sumassets {
    qui reg `Y' c.l1lnaid_wb##i.treat_old ///
        $rddaidcontrol fdist_viloldb seg_viloldb2dum* ///
        if year==2005 & dist_viloldb<50 & province!="bamyan" & province!="paktya" & province!="kabul", ///
        cluster(districtcode)
    sum `Y' if e(sample)
    cap drop `Y'_sd
    gen double `Y'_sd = `Y' / r(sd)
}


matrix MC = J(6,2,.)
matrix coln MC = b se
matrix rown MC = "(log) Night light" "Wheat consumption" "(log) Food expenditure" "Dietary diversity" "Food security" "Sum of assets"

matrix MT = J(6,2,.)
matrix coln MT = b se
matrix rown MT = "(log) Night light" "Wheat consumption" "(log) Food expenditure" "Dietary diversity" "Food security" "Sum of assets"

local n 1

* --- NIGHTLIGHT (district level) ---
preserve
    keep if year==2005 & dist_viloldb<45 & province!="bamyan" & province!="paktya" & province!="kabul"
    
    * Collapse to district level
    collapse (mean) log_nightlight l1lnaid_wb_sd treat_old $rddaidcontrol fdist_viloldb seg_viloldb2dum*, ///
        by(districtcode)
    
    * Standardize outcome at district level
    quietly sum log_nightlight
    gen double log_nightlight_sd = log_nightlight / r(sd)
    
    qui reg log_nightlight_sd c.l1lnaid_wb_sd##i.treat_old ///
        $rddaidcontrol fdist_viloldb seg_viloldb2dum*, ///
        cluster(districtcode)
    
    * Control group (treat_old = 0): ME = beta_aid
    matrix MC[`n',1] = _b[l1lnaid_wb_sd]
    matrix MC[`n',2] = _se[l1lnaid_wb_sd]
    
    * Treatment group (treat_old = 1): ME = beta_aid + interaction
    lincom l1lnaid_wb_sd + 1.treat_old#c.l1lnaid_wb_sd
    matrix MT[`n',1] = r(estimate)
    matrix MT[`n',2] = r(se)
restore

local n = `n' + 1

* --- Household-level outcomes ---
foreach Y of varlist amount_wheat_flour_sd lnhhexp_totp2005_sd dietarydiversity_sd foodsecure_sd sumassets_sd {
    qui reg `Y' c.l1lnaid_wb_sd##i.treat_old ///
        fdist_viloldb seg_viloldb2dum* $rddaidcontrol ///
        if year==2005 & dist_viloldb<45 & province!="bamyan" & province!="paktya" & province!="kabul", ///
        cluster(districtcode)
    
    * Control group
    matrix MC[`n',1] = _b[l1lnaid_wb_sd]
    matrix MC[`n',2] = _se[l1lnaid_wb_sd]
    
    * Treatment group
    lincom l1lnaid_wb_sd + 1.treat_old#c.l1lnaid_wb_sd
    matrix MT[`n',1] = r(estimate)
    matrix MT[`n',2] = r(se)
    
    local n = `n' + 1
}


coefplot ///
    (matrix(MC[,1]), se(MC[,2]) label("Control") ///
        mcolor(gs4) ciopts(lcolor(gs4 gs4)) levels(99 95)) ///
    (matrix(MT[,1]), se(MT[,2]) label("ISAF Treatment") ///
        msymbol(D) mcolor(gs10) ciopts(lcolor(gs10 gs10)) levels(99 95)) ///
    , ///
    xline(0, lcolor(cranberry) lwidth(thin) lpattern(solid)) ///
	xtitle("Effect on standardized outcome", margin(medium)) ///
    coeflabels( ///
		1 = "(log) Night light" ///
		2 = "Wheat consumption" ///
        3 = "(log) Food expenditure" ///
        4 = "Dietary diversity" ///
        5 = "Food security" ///
        6 = "Sum of assets" ///
        , wrap(40) truncate(40) notick labsize(medium) labgap(3) ///
    ) ///
    graphregion(color(white)) ///
    legend(position(7) rows(2) region(lstyle(black)))
	
graph save "graphs\aidcoefplot.gph", replace
graph export "graphs\aidcoefplot.pdf", replace		
		
			

/****************************************************************************************
* FIGURE 5 - Panel: Heterogeneity in aid effectiveness, 2005-2010
****************************************************************************************/

************ Panel A: Interaction effect (Any base x WB AID) *******************

cap grstyle init
cap grstyle color background white
cap mkdir figures

local outcomes ///
    log_nightlight amount_wheat_flour lnhhexp_totp2005 ///
    dietarydiversity foodsecure sumassets

local samp "year>2003 & year<2011"

local fe_hh "id i.year"
local fe_d  "districtcode i.year"

local K : word count `outcomes'
matrix I_cur = J(`K',2,.)
matrix I_lag = J(`K',2,.)

local coeflabs ""
local no = 1

foreach y of local outcomes {

    * Y-axis label
    local lab : variable label `y'
    if "`lab'"=="" local lab "`y'"
    local coeflabs `"`coeflabs' r`no'="`lab'""'

    *==========================
    * DISTRICT-LEVEL outcome
    *==========================
    if inlist("`y'", "log_nightlight") {

        preserve
            keep if `samp'

            * Collapse to district-year with all vars needed for BOTH regressions
            collapse (mean) ///
                `y' activebase l1activebase ///
                lnaid_wb l1lnaid_wb ///
                l1lnallsigacts l1vhi ///
                , by(districtcode year)

            * Standardize outcome at district-year level
            quietly summarize `y'
            gen z_`y' = (`y' - r(mean)) / r(sd)

            * --- Current: Any base × Aid (current) ---
            quietly reghdfe z_`y' i.activebase##c.lnaid_wb, ///
                absorb(`fe_d') cluster(districtcode)
            matrix I_cur[`no',1] = _b[1.activebase#c.lnaid_wb]
            matrix I_cur[`no',2] = _se[1.activebase#c.lnaid_wb]

            * --- Lagged: Any base × Aid (lagged) ---
            quietly reghdfe z_`y' i.l1activebase##c.l1lnaid_wb, ///
                absorb(`fe_d') cluster(districtcode)
            matrix I_lag[`no',1] = _b[1.l1activebase#c.l1lnaid_wb]
            matrix I_lag[`no',2] = _se[1.l1activebase#c.l1lnaid_wb]

        restore
    }

    *==========================
    * HOUSEHOLD-LEVEL outcomes
    *==========================
    else {

        * Ensure standardized outcome exists
        cap confirm variable z_`y'
        if _rc {
            di as error "Missing z_`y'. Please standardize `y' first (create z_`y')."
            exit 198
        }

        * --- Current: Any base × Aid (current) ---
        quietly reghdfe z_`y' i.activebase##c.lnaid_wb ///
            [pweight=hhweight] if `samp', absorb(`fe_hh') cluster(districtcode)
        matrix I_cur[`no',1] = _b[1.activebase#c.lnaid_wb]
        matrix I_cur[`no',2] = _se[1.activebase#c.lnaid_wb]

        * --- Lagged: Any base × Aid (lagged) ---
        quietly reghdfe z_`y' i.l1activebase##c.l1lnaid_wb ///
            [pweight=hhweight] if `samp', absorb(`fe_hh') cluster(districtcode)
        matrix I_lag[`no',1] = _b[1.l1activebase#c.l1lnaid_wb]
        matrix I_lag[`no',2] = _se[1.l1activebase#c.l1lnaid_wb]
    }

    local no = `no' + 1
}


coefplot ///
    (matrix(I_cur[,1]), se(I_cur[,2]) label("Contemporaneous") pstyle(p2) msymbol(O) mcolor(black) ciopts(color(black black)) level(99 95)) ///
    (matrix(I_lag[,1]), se(I_lag[,2]) label("Lagged")  pstyle(p5) msymbol(D) mcolor(gs8) ciopts(color(gs8 gs8)) level(99 95)) ///
    , ///
    xtitle("Change in aid effect", margin(medium)) ///
    xline(0, lwidth(thin) lpattern(solid)) ///
    xlabel(-0.8(0.2)0.6, nogrid labgap(1)) ///
    coeflabels(`coeflabs') ///
    legend(position(7) rows(1) region(lstyle(black))) ///
    scale(0.8) ///
    graphregion(color(white))

graph export "graphs\lnaid_wb_anybase.pdf", replace
graph save   "graphs\lnaid_wb_anybase.gph", replace


*********** Panel B: Interaction effect (Any base x PRT projects)) *************

cap grstyle init
cap grstyle color background white
cap mkdir figures

local outcomes ///
    log_nightlight amount_wheat_flour lnhhexp_totp2005 ///
    dietarydiversity foodsecure sumassets

local samp "year>2003 & year<2011"

local fe_hh "id i.year"
local fe_d  "districtcode i.year"

*-----------------------------
* 1) SD of PRT projects (current), in estimation sample
*-----------------------------
quietly summarize no_prtproject if `samp'
local sd_prtproj = r(sd)
if missing(`sd_prtproj') | `sd_prtproj'==0 {
    di as error "SD(no_prtproject) is missing/zero in the estimation sample. Cannot SD-scale."
    exit 459
}
di as txt "SD(no_prtproject) in sample = " as res `sd_prtproj'

local K : word count `outcomes'
matrix I_cur = J(`K',2,.)
matrix I_lag = J(`K',2,.)

local coeflabs ""
local no = 1

foreach y of local outcomes {

    local lab : variable label `y'
    if "`lab'"=="" local lab "`y'"
    local coeflabs `"`coeflabs' r`no'="`lab'""'

    *==========================
    * DISTRICT-LEVEL outcome
    *==========================
    if "`y'"=="log_nightlight" {

        preserve
            keep if `samp'

            * Collapse with all vars needed
            collapse (mean) ///
                `y' activebase l1activebase ///
                no_prtproject l1no_prtproject ///
                l1lnallsigacts l1vhi log_l1nightlight loan_i ///
                , by(districtcode year)

            quietly summarize `y'
            gen z_`y' = (`y' - r(mean)) / r(sd)

            * Current interaction
            quietly reghdfe z_`y' i.activebase##c.no_prtproject, ///
                absorb(`fe_d') cluster(districtcode)
            matrix I_cur[`no',1] = _b[1.activebase#c.no_prtproject]  * `sd_prtproj'
            matrix I_cur[`no',2] = _se[1.activebase#c.no_prtproject] * `sd_prtproj'

            * Lagged interaction (scaled by SD of current no_prtproject)
            quietly reghdfe z_`y' i.l1activebase##c.l1no_prtproject, ///
                absorb(`fe_d') cluster(districtcode)
            matrix I_lag[`no',1] = _b[1.l1activebase#c.l1no_prtproject]  * `sd_prtproj'
            matrix I_lag[`no',2] = _se[1.l1activebase#c.l1no_prtproject] * `sd_prtproj'

        restore
    }

    *==========================
    * HOUSEHOLD-LEVEL outcomes
    *==========================
    else {

        cap confirm variable z_`y'
        if _rc {
            di as error "Missing z_`y' -> standardize `y' first."
            exit 198
        }

        quietly reghdfe z_`y' i.activebase##c.no_prtproject  ///
            [pweight=hhweight] if `samp', absorb(`fe_hh') cluster(districtcode)
        matrix I_cur[`no',1] = _b[1.activebase#c.no_prtproject]  * `sd_prtproj'
        matrix I_cur[`no',2] = _se[1.activebase#c.no_prtproject] * `sd_prtproj'

        quietly reghdfe z_`y' i.l1activebase##c.l1no_prtproject ///
            [pweight=hhweight] if `samp', absorb(`fe_hh') cluster(districtcode)
        matrix I_lag[`no',1] = _b[1.l1activebase#c.l1no_prtproject]  * `sd_prtproj'
        matrix I_lag[`no',2] = _se[1.l1activebase#c.l1no_prtproject] * `sd_prtproj'
    }

    local no = `no' + 1
}


coefplot ///
    (matrix(I_cur[,1]), se(I_cur[,2]) label("Contemporaneous") pstyle(p2) msymbol(O) mcolor(black) ciopts(color(black black)) level(99 95)) ///
    (matrix(I_lag[,1]), se(I_lag[,2]) label("Lagged")  pstyle(p5) msymbol(D) mcolor(gs8) ciopts(color(gs8 gs8)) level(99 95)) ///
    , ///
    xtitle("Change in PRT-project effect", margin(medium)) ///
    xline(0, lwidth(thin) lpattern(solid)) ///
	xlabel(-0.8(0.2)0.6, nogrid labgap(1)) ///
    coeflabels(`coeflabs') ///
    legend(position(7) rows(1) region(lstyle(black))) ///
    scale(0.8) ///
    graphregion(color(white))

graph export "graphs\no_prtproject_anybase.pdf", replace
graph save   "graphs\no_prtproject_anybase.gph", replace


		
********************************************************************************
* TABLE 1 - GRD: Different outcomes (NRVA), 2005
********************************************************************************

cap drop indivbw
gen indivbw = .

foreach X of varlist copehelp copesocial {

    rdbwselect `X' fdist_viloldb anyshock $rddcontrols dum_segment200km* ///
        if year==2005 & province!="bamyan" & province!="paktya" & province!="kabul", ///
        c(0) p(1) kernel(uniform) vce(cluster id)

    replace indivbw = e(h_mserd)

    reghdfe `X' treat_old fdist_viloldb anyshock $rddcontrols [pweight=hhweight] ///
        if year==2005 & dist_viloldb<indivbw & province!="bamyan" & province!="paktya" & province!="kabul", ///
        absorb(segment200km) cluster(districtcode)
    estadd local seg "Yes"
    estadd local indvw = round(indivbw, .1)
    estadd local control "Yes"
    estadd local restr "Yes"
    est sto rdd1`X'

    reghdfe `X' i.treat_old##c.llnwheat_rw_shock fdist_viloldb anyshock $rddcontrols [pweight=hhweight] ///
        if year==2005 & dist_viloldb<indivbw & province!="bamyan" & province!="paktya" & province!="kabul", ///
        absorb(segment200km) cluster(districtcode)
    estadd local seg "Yes"
    estadd local indvw = round(indivbw, .1)
    estadd local control "Yes"
    estadd local restr "Yes"
    est sto rdd3`X'

    reghdfe `X' i.treat_old##c.l1lnallsigacts fdist_viloldb anyshock $rddcontrols [pweight=hhweight] ///
        if year==2005 & dist_viloldb<indivbw & province!="bamyan" & province!="paktya" & province!="kabul", ///
        absorb(segment200km) cluster(districtcode)
    estadd local seg "Yes"
    estadd local indvw = round(indivbw, .1)
    estadd local control "Yes"
    estadd local restr "Yes"
    est sto rdd4`X'

    reghdfe `X' i.treat_old##i.shock fdist_viloldb anyshock $rddcontrols [pweight=hhweight] ///
        if year==2005 & dist_viloldb<indivbw & province!="bamyan" & province!="paktya" & province!="kabul", ///
        absorb(segment200km) cluster(districtcode)
    estadd local seg "Yes"
    estadd local indvw = round(indivbw, .1)
    estadd local control "Yes"
    estadd local restr "Yes"
    est sto rdd5`X'
}


foreach X of varlist helpsocial shura_member cdc_member {

    rdbwselect `X' fdist_viloldb $rddcontrols_short dum_segment200km* ///
        if year==2005 & province!="bamyan" & province!="paktya" & province!="kabul", ///
        c(0) p(1) kernel(uniform) vce(cluster id)

    replace indivbw = e(h_mserd)

    reghdfe `X' treat_old fdist_viloldb $rddcontrols_short [pweight=hhweight] ///
        if year==2005 & dist_viloldb<indivbw & province!="bamyan" & province!="paktya" & province!="kabul", ///
        absorb(segment200km) cluster(districtcode)
    estadd local seg "Yes"
    estadd local indvw = round(indivbw, .1)
    estadd local control "Yes"
    estadd local restr "Yes"
    est sto rdd1`X'

    reghdfe `X' i.treat_old##c.llnwheat_rw_shock fdist_viloldb $rddcontrols_short [pweight=hhweight] ///
        if year==2005 & dist_viloldb<indivbw & province!="bamyan" & province!="paktya" & province!="kabul", ///
        absorb(segment200km) cluster(districtcode)
    estadd local seg "Yes"
    estadd local indvw = round(indivbw, .1)
    estadd local control "Yes"
    estadd local restr "Yes"
    est sto rdd3`X'

    reghdfe `X' i.treat_old##c.l1lnallsigacts fdist_viloldb $rddcontrols_short [pweight=hhweight] ///
        if year==2005 & dist_viloldb<indivbw & province!="bamyan" & province!="paktya" & province!="kabul", ///
        absorb(segment200km) cluster(districtcode)
    estadd local seg "Yes"
    estadd local indvw = round(indivbw, .1)
    estadd local control "Yes"
    estadd local restr "Yes"
    est sto rdd4`X'

    reghdfe `X' i.treat_old##i.shock fdist_viloldb $rddcontrols_short [pweight=hhweight] ///
        if year==2005 & dist_viloldb<indivbw & province!="bamyan" & province!="paktya" & province!="kabul", ///
        absorb(segment200km) cluster(districtcode)
    estadd local seg "Yes"
    estadd local indvw = round(indivbw, .1)
    estadd local control "Yes"
    estadd local restr "Yes"
    est sto rdd5`X'
}

			esttab rdd1copehelp rdd1copesocial rdd1helpsocial rdd1shura_member rdd1cdc_member using "tables\RDD_A_alloutcomes.tex", fragment booktabs se(3) b(3) star starlevels(* 0.10 ** 0.05 *** 0.01) replace  ///
					mlabel(none) style(fixed)  keep(treat_old) order(treat_old) ///
					varlabel( treat_old "ISAF" ) ///
					stats(r2_a N  N_clust indvw , fmt(3 0 0 3  ) label("Adj. R$^2$" "\hline  Observations"  "Districts" "Bandwidth"  )) nonumbers nogaps noabbr nolines nomtitle nonote 
					
									
********************************************************************************
* TABLE 2 - Panel: Community help (NRVA), 2005-2008
********************************************************************************

	cap est drop *
	*Panel A
		foreach X of varlist largebase_inbr activebase  no_securityprt {
			replace contest=l1lnallsigacts
			replace isaf=`X'
						
			reghdfe copehelp isaf contest  $short anyshock [pweight=hhweight] if year>2003  & year<2011  , absorb(id i.year ) cluster(districtcode) 
			estadd local controls "Yes"
			estadd local lfe "Yes"
			estadd local lag "No"
			est store b1`X'
			
			replace isaf=l1`X'
			*replace contest=l2lnallsigacts			
			reghdfe copehelp isaf contest  $short  anyshock [pweight=hhweight] if year>2003 & year<2011, absorb(id i.year ) cluster(districtcode) 
			estadd local controls "Yes"
			estadd local lfe "Yes"
			estadd local lag "Yes"
			est store b2`X'
			
			}
		
		
		esttab b1largebase_inbr b2largebase_inbr b1activebase b2activebase b1no_securityprt b2no_securityprt using "tables\PanelAMain.tex", fragment booktabs se(3) b(3) star starlevels(* 0.10 ** 0.05 *** 0.01) replace  ///
			mlabel(none) style(fixed) keep(isaf ) order(isaf shockisaf ) ///
			varlabel( isaf "ISAF" ) ///
			stats(r2_a N, fmt(3 0) label("Adj. R$^2$" "Observations")) nonumbers nogaps noabbr nomtitle nolines nonote
		
			
	*Panel B
	cap est drop *
		
		foreach X of varlist largebase_inbr activebase  no_securityprt {
			replace contest=l1lnallsigacts
			replace isaf=`X'
						
			reghdfe copehelp i.isaf##c.llnwheat_rw_shock anyshock   $short   [pweight=hhweight] if year>2003  & year<2011, absorb(id i.year ) cluster(districtcode) 
			estadd local controls "Yes"
			estadd local lfe "Yes"
			estadd local lag "No"
			est store b1`X'

			replace isaf=l1`X'
						
			reghdfe copehelp i.isaf##c.llnwheat_rw_shock  anyshock  $short   [pweight=hhweight] if year>2003  & year<2011 , absorb(id i.year ) cluster(districtcode) 
			estadd local controls "Yes"
			estadd local lfe "Yes"
			estadd local lag "Yes"
			est store b2`X'
			
			
			}
		
		esttab b1largebase_inbr b2largebase_inbr b1activebase b2activebase b1no_securityprt b2no_securityprt using "tables\PanelBMain.tex", fragment booktabs se(3) b(3) star starlevels(* 0.10 ** 0.05 *** 0.01) replace  ///
			mlabel(none) style(fixed) keep(1.isaf 1.isaf#c.llnwheat_rw_shock llnwheat_rw_shock) order(1.isaf 1.isaf#c.llnwheat_rw_shock llnwheat_rw_shock) ///
			varlabel( 1.isaf "ISAF" 1.isaf#c.llnwheat_rw_shock "ISAF*Shock" llnwheat_rw_shock "Shock") ///
								stats(r2_a N lag, fmt(3 0 0) label("Adj. R$^2$" "Observations" "\hline  Lagged")) nonumbers nogaps noabbr nomtitle nolines nonote
	
	
	
	
********************************************************************************
* TABLE 7 - GRD: Mechanisms (NRVA), 45km bandwidth, 2005
********************************************************************************
				
				cap est drop *

				foreach X of varlist copemilstate loan govemployment agrincome eradication {
				
				reghdfe `X' treat_old fdist_viloldb [pweight=hhweight] if year==2005 & dist_viloldb<45 &  province!="bamyan" & province!="paktya" & province!="kabul", absorb(segment200) cluster(districtcode)
					estadd local seg "Yes"
					estadd local control "No"
					estadd local restr "Yes"
					est store chan1`X'
				}
				
				foreach X of varlist f_enrolled m_enrolled sh_kids_enroll sh_f_enroll sh_f_m_enroll {
				
				reghdfe `X' treat_old fdist_viloldb [pweight=hhweight] if year==2005 & dist_viloldb<45 &  province!="bamyan" & province!="paktya" & province!="kabul", absorb(segment200) cluster(districtcode)
					estadd local seg "Yes"
					estadd local control "No"
					estadd local restr "Yes"
					est store chan2`X'
				}
				
				
				*With controls: govemployment, eradication turns insignificant
				
				foreach X of varlist amount_wheat_flour hhexp_totp2005 dietarydiversity foodinsecure sumassets  {
				reghdfe `X' treat_old fdist_viloldb [pweight=hhweight] if year==2005 & dist_viloldb<45 &  province!="bamyan" & province!="paktya" & province!="kabul", absorb(segment200) cluster(districtcode)
					estadd local seg "Yes"
					estadd local control "No"
					estadd local restr "Yes"
					est store chan3`X' 
				}
				
				* anyaid_comm anyaid_disb lnaid_afg lnaid_wb no_project
				foreach X of varlist cashforwork any_cdc lnaid_wb log_nightlight economicimprove {
				reghdfe `X' treat_old fdist_viloldb [pweight=hhweight] if year==2005 & dist_viloldb<45 &  province!="bamyan" & province!="paktya" & province!="kabul", absorb(segment200) cluster(districtcode)
					estadd local seg "Yes"
					estadd local control "No"
					estadd local restr "Yes"
					est store chan4`X' 
				}

				
				*counterinsurgency
				foreach X of varlist dis_insecurity lnallsigacts lnbest_est lndeaths_civilians lndeaths_civilianssb {
				reghdfe `X' treat_old fdist_viloldb [pweight=hhweight] if year==2005 & dist_viloldb<45 &  province!="bamyan" & province!="paktya" & province!="kabul", absorb(segment200) cluster(districtcode)
					estadd local seg "Yes"
					estadd local control "No"
					estadd local restr "Yes"
					est store chan5`X' 
				}
				
									
				esttab chan1*  using "tables\channel1.tex", fragment booktabs se(3) b(3) star starlevels(* 0.10 ** 0.05 *** 0.01) replace  ///
						mlabel(none) style(fixed)  keep(treat_old) order(treat_old) ///
						varlabel( treat_old "ISAF" ) ///
						stats(r2_a N, fmt(3 0) label("Adj. R$^2$" "Observations")) nonumbers nogaps noabbr nolines nomtitle nonote postfoot( \hline)
				
				esttab chan2*  using "tables\channel2.tex", fragment booktabs se(3) b(3) star starlevels(* 0.10 ** 0.05 *** 0.01) replace  ///
						mlabel(none) style(fixed)  keep(treat_old) order(treat_old) ///
						varlabel( treat_old "ISAF" ) ///
						stats(r2_a N, fmt(3 0) label("Adj. R$^2$" "Observations")) nonumbers nogaps noabbr nolines nomtitle nonote postfoot( \hline)
							
				esttab chan3*  using "tables\channel3.tex", fragment booktabs se(3) b(3) star starlevels(* 0.10 ** 0.05 *** 0.01) replace  ///
						mlabel(none) style(fixed)  keep(treat_old) order(treat_old) ///
						varlabel( treat_old "ISAF" ) ///
						stats(r2_a N, fmt(3 0) label("Adj. R$^2$" "Observations")) nonumbers nogaps noabbr nolines nomtitle nonote postfoot( \hline)
				
				esttab chan4*  using "tables\channel4.tex", fragment booktabs se(3) b(3) star starlevels(* 0.10 ** 0.05 *** 0.01) replace  ///
						mlabel(none) style(fixed)  keep(treat_old) order(treat_old) ///
						varlabel( treat_old "ISAF" ) ///
						stats(r2_a N seg control restr, fmt(3 0 0 0 0 ) label("Adj. R$^2$" "Observations" "\hline  200km segments" "Control variables" "Restricted sample")) nonumbers nogaps noabbr nolines nomtitle nonote 
	
				esttab chan5*  using "tables\channel5.tex", fragment booktabs se(3) b(3) star starlevels(* 0.10 ** 0.05 *** 0.01) replace  ///
						mlabel(none) style(fixed)  keep(treat_old) order(treat_old) ///
						varlabel( treat_old "ISAF" ) ///
						stats(r2_a N, fmt(3 0) label("Adj. R$^2$" "Observations")) nonumbers nogaps noabbr nolines nomtitle nonote postfoot( \hline)
	
	
						est drop * 


